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I 



INTRODUCTION 



Computational fluid dynamics has matured significantly 
within the past decade because of the development of 
increased computational capabilities and powerful 
computational techniques. Current problems being addressed 
include predicting flow separation over airfoils and 
post-stall flight characteristics. These areas are of 
interest because studies [Ref. 1] indicate increased lift 
and thus sustained flight are attainable when an airfoil is 
dynamically stalled, that is, its angle of attack is pitched 
to a post stall angle of attack rather than being initially 
placed at that high lift. 

Computational methods utilizing the full Navier-Stokes 
(N-S) equations are capable of addressing these issues, as 
are methods that include approximations to the Navier-Stokes 
equations. One method, the Interactive Boundary Layer (IBL) 
technique, developed by Tuncer Cebeci at Douglas Aircraft 
Company and at the California State University [Ref. 2], 
divides the flow over an airfoil into a viscous inner 
boundary layer and an inviscid outer layer. The 
characteristics of the inner flow are obtained from a 
numerical solution of Prandtl • s boundary layer equation and 
the outer flow's characteristics are determined from Hess 
and Smith's panel method, and Fourier analysis and conformal 
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mapping. The inner and outer layers are then redetermined 
by an interaction model that iterates between the two 
regions and marches downstream until the flow conditions 
have been satisfied at the boundary for both regions. The 
Cebeci IBL code uses Michel's criterion to predict 
transition from laminar to turbulent flow or transition may 
be prescribed. An algebraic (Cebeci-Smith) turbulence model 
is used. 

A full Navier-Stokes code developed by N.L. Sankar and 
his associates at the Georgia Institute of Technology [Ref. 
3] uses an implicit finite-difference procedure to solve the 
2-D Reynolds-averaged compressible Navier-Stokes equations 
in strong conservative form. The time-marching algorithm 
used is an Alternating Direction Implicit (ADI) procedure 
developed by Beam and Warming [Ref. 4] and implemented by 
Steger [Ref. 5] . The Sankar N-S code uses a body-fitted 
C-grid system and an algebraic (Baldwin-Lomax) turbulence 
model . 

The Cebeci IBL and the Sankar N-S codes are designed for 
different purposes. A low Reynolds number flow over an 
airfoil tends to be laminar until separation. The flow then 
transitions to turbulent flow and reattaches as turbulent 
flow. The Cebeci IBL code models this separation bubble if 
transition is specified within the separation bubble. 
Velocity profiles and skin coefficients are extremely 
important in analyzing these low Reynolds flows. Cebeci has 
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developed codes for compressible oscillating airfoils, 
however, this Cebeci IBL code was developed for 
incompressible steady state flow only and thus does not 
predict the effects of unsteady flow nor compressibility. 

The Sankar N-S code was developed to address dynamic 
stall and its implication of increased lift. Therefore, the 
values of interest to date have been coefficients of 
pressure, lift, moment, and the effect of hysteresis on 
these values. However, the Sankar N-S code assumes the flow 
is fully turbulent, and therefore does not account for 
transition from laminar to turbulent flow. 

Neither dynamic stall nor transition within a separation 
bubble are easily quantified experimentally. Transition is 
a boundary layer phenomenon and the velocity profiles and 
skin frictions within the boundary layer must be measured to 
assure correct interpretation of the flow under 
investigation; surface pressures are not sufficient to 
accurately locate flow separation and reattachment [Ref. 6]. 
This is a time consuming, expensive process prone to error. 
Experimental methods include laser anemometry and hot wire 
probes [Ref. 7]. Disturbance of the boundary layer flow due 
to probes is undesirable and hot wire probes are normally 
unable to determine flow direction; therefore laser 
anemometry, although expensive and tedious, is increasingly 
being used. 
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Dynamic stall is difficult to characterize due to the 
transitory nature of the phenomenon. Experimental 
techniques and apparatus include pressure transducers, hot 
wire probes, and laser doppler velocimetry [Ref. 8]. Flow 
visualization is also a very effective tool for studying 
both dynamic stall and separation bubbles [Refs. 9,10]. 

A good review of the current state-of-the-art 
computational and experimental aspects of aerodynamic flows 
is given in the proceedings of three symposia on this topic 
edited by T. Cebeci [Ref. 11]. 
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II 



OBJECTIVES 



The intent of this study is two-fold: to become 
familiar with computational fluid dynamic methods and to 
evaluate two codes to determine their range of 
applicability. 

Computational fluid dynamics consists of various 
mathematical methods and implementation schemes. A 
significant portion of analysis inherent in computational 
codes is empirical; therefore, the assumptions used strongly 
influence the results. It is important, when attempting to 
choose a computational code for a specific purpose, to be 
familiar with the significance of the analytical methods, 
assumptions made, and empirical models. Each code is 
different in these respects and must be analyzed 
individually and in detail to assure reliable, accurate 
results, especially when extending the flow regime or 
airfoil to conditions whose features are unknown. 

The two codes chosen, the Cebeci IBL code and the Sankar 
N-S code, are a good representation of two powerful, 
accurate methods that differ widely in computational 
approach. The Cebeci IBL code has been extensively tested 
in a variety of steady state conditions [Ref. 12] and the 
Sankar N-S code has compared well to experimental data for a 
pitching airfoil [Ref. 13]. However, the lack of steady 
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state boundary layer data for the Sankar N-S code indicated 
that a more in-depth analysis of the applicability of the 
code was required. The Cebeci IBL code was used as the 

reference for the Sankar code. 

The analysis included the following: 

1. Assess Ci for Reynolds numbers of 1.5 and 6 million 

2. Assess Cp and Cf for a range of Reynolds numbers (1-15 
million at 0 degrees angle of attack) and angles of 
attack ( 0 , 2 , 4 , and 6 degrees for 1 . 5M Re number and 
0, 4, 8 , and 12 degrees for 6M Re number) 

3. Assess velocity profiles for Reynolds numbers of 1, 6, 
and 15 million 

4. Assess the influence of dissipation factors and 
grid size on the results of the Sankar N-S code. 
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III. MATHEMATICAL FORMULATION 



A. GOVERNING EQUATIONS 

Flow over an airfoil can be described by the velocity 

A A A A 

q = ui + vj + wk, the pressure, the density, and the 
temperature. These six variables (u, v, w, p, p, and T) are 
fully described by the continuity equation, the equation of 
state, p = pRT; the energy equation, 6Q - 6W = 6E; and the 
three equations of motion. [Ref. 14] 

The continuity equation states mass is conserved; 
i.e., the flux of mass through a cube per time is equal to 
the time rate of change of mass. This is shown in Figure 
3.1 for the flow through the faces perpendicular to the x 
axis. Mathematically this is expressed as 

[ 3 ^ U - Ax] AyAz - [ 3 ^ v) AylAzAx - [ - ^ Az] AxAy 

= (pAxAyAz) (3.1) 

Since the control volume is fixed, AxZyAz is 
independent of time; therefore 

3a + iiM + iigu o , (3.2, 

or 
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Net flux through face perpendicular to x axis 

_ _ r 9 ( pu )_ Ax j Ay A z 
L 3x 

Source: [Ref. 14:p. 106] 

Figure 3 . 1 The Flux of Mass Through a Cube 



Ax] AyAz 
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3 +: + A • pq - 0 . 



(3.3) 



Alternately, 



Dp ,3u 3v 3w. _ n 
Dt + p fe + 3y + 3i- } " ° 



(3.4) 



or 



§£ + p(A-q) = 0 



(3.5) 



The three equations of motion, one for each axis of the 
Cartesian coordinate system, are described by Newton's 
second law, AF = A (ma) . The summation of the x components 
of surface forces on the element shown in Figure 3.2 is 



AF = Ama = (pAxAyAz) 

X X 



Du 

Dt 



(3.6) 



with 



3a 

AF x = - a x AyAz + ■ (a x + Ax) AyAz 

9 V 

- Tyx AxAz + (iy X + Ay) AxAz 

3x 

zx 

- •'zx toiy + (T zx + — 4z)AxAy • 
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STRESS-STRAIN RELATIONS 




Figure 3.2 X-components of Surface Forces on an Element 
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Dividing by the volume of the element yields 



8a 8 t 8t 

x . yx . zx 



8x 8y 



8z 



Du 
P Dt 



(3.8) 



Similarly, for the y and z directions 



8a 8 t 8t _ 

_y + x y +■... z y_ = D dv 

8y 8x 8z p Dt 



(3.9) 



and 



8a 8t 8t _ 

z xz yz Dw 

8z 8x 8y P Dt 



(3.10) 



For Newtonian fluids with a single viscosity coefficient, 
the normal and tangential shear stresses are as follows 
[Ref. 15]: 



° x = -p + 2lJ li - (V '5> 



(3.11a) 



a = - 



p + 2y fy " (V ’^ 



(3.11b) 



a_ = - 



P + 2 V §7 - |y (V-q) 



(3.11c) 



T = T 
yx xy 



,9v , 3u. 

u 1 35 + 5y > 



(3. lid) 



T = T 
yz zy 



,8w 8Vv 

w( 57 + 3i> 



(3. lie) 



,8u , 8w. 

T = T = ]J (tt— + trr?) 

zx xy 8z 8x 



(3. Ilf) 
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Substituting these stress definitions into the equations of 
motion, and assuming a constant viscosity corresponding to 
the mean temperature of the fluid ultimately yields the 
Navier-Stokes equations: 



3x 



3y 



3z 



y [ 



,2 




,.2 




.2 


3 u 




3 u 




3 u 


„ 2 


+ 


■v 2 


+ 


* 2 


3x 




3y 




3z 


*2 
3 v 




~ 2 
3 v 




3 2 v 


2 


+ 


„ 2 


+ 


„ 2 


3x 




3y 




3z 


~ 2 
3 w 




3 2 w 




3 2 w 


„ 2 


+ 


2 


+ 


„ 2 


3x 




3y 




3z 



3 3x 



] + $ |y[7-q] 



3 3z 



= P 



= P 



Du 

Dt 



Dv 

Dt 



„ Dw 
p Dt 



(3.12a) 



(3.12b) 



(3.12c) 



or in vector format 



-Vp + yV 2 q + yV ( V • q ) 



p |3. + p (q • V ) q 



(3.13) 
[Ref. 14] 



B. REYNOLDS STRESSES 

The Navier-Stokes equations are valid for laminar and 
turbulent flow. However, the complexity of turbulence has 
made it impossible to relate the motion of the fluid to the 
boundary conditions and obtain an exact solution. 
Therefore, the turbulence must currently be modeled. 0. 
Reynolds divided the turbulent flow into a mean motion and 
fluctuating, or eddying, motion as follows: 
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u 



u + u' 



(3.14a) 



V = V + V ' 



(3.14b) 



w = w + w ' 



(3.14c) 



P = P + P' 



( 3 . 14d) 



P = P + P' 



(3 . 14e) 



T = T + T' 



(3 . 14f ) 



where the barred terms are the time-average of the component 
and the slashed terms are the fluctuations. By definition, 
the time averages of all quantities describing the 
fluctuations are equal to zero: 

u' =0, v' = 0, w ' = 0, p ' = 0, p"' = 0, T ' = 0 

Rules for operating on mean time-averages are given 
below. F and g are dependent variables, and s is the 
independent variable x, y, z, or t. 



f = f 



f + g = TT + g 



f * g = f * g. 



[Ref. 16] 
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The stresses caused by the fluctuations can be 
determined using the momentum theorem. Consider an area dA 
with dA * pu * dt being the mass of incompressible fluid 
passing through the element in time dt. Thus, the flux of 
momentum in the x direction is 

dJ x = dA * pu 2 *dt ; (3.15a) 



Correspondingly , 



dJy = dA * p uv * dt 



(3.15b) 



and 



dJ z = dA 


* p uw * dt. 


(3.15c) 


Calculating the time 


averages for the 


fluxes of momentum 


unit time yields: 


— ^ 




dJ x = dA 


P 


(3.16a) 








s 

II 

>1 


p uv 


(3.16b) 


c[J z = dA 


p uw. 


(3.16c) 
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Utilizing the definition of turbulent flow and the previous 
rules yields 



dJ = dA • p (u^ + u'2) 



(3.17a) 



dJy = dA • p(u*v + u'v 1 ) 



(3.17b) 



dJ = dA • p (u*w + u'w') 



(3.17c) 



Dividing these rates of change of momentum by area dA, we 
obtain stresses. The equal and opposite stresses exerted on 
the area by the surroundings are a normal stress, -(u 2 + 
u~' 2 ), and two shearing stresses, -(uv + u'v' ) and -(uw + 
u'w') . Thus, the superposition of fluctuations on the mean 
motion gives rise to three additional stresses 
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-pu'v' , 



t' 



xz 



-pu'w' 



(3.18) 



The total stress tensor due to the turbulent velocity 
components of the flow are 
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The presence of fluctuations presents itself as an 
apparent increase in stresses (viscosity) . These 
additional stresses over the mean or laminar stresses are 
termed apparent stresses or Reynolds stresses. 

B. TURBULENCE MODELING 

The presence of Reynolds stresses in turbulent flow 
introduces additional unknowns in the Navier-Stokes 
equations. Therefore, the Navier-Stokes equations, 
continuity, the perfect gas law, and the energy equation are 
no longer sufficient to completely define a solution. This 
is known as the closure problem and is usually resolved by 
turbulence modeling. [Ref. 17] 

A common method used is to relate the turbulent stress 
to the mean flow properties through empirically based 
algebraic formulas. An eddy viscosity, is defined in 
the same form as the laminar viscosity. Previous models 
related surface boundary conditions to points in the fluid 
away from the boundaries through wall functions. This 
avoided modeling the direct influence of the eddy viscosity; 
however, it is only applicable in regions where the Reynolds 
number is high enough for viscous effects to be unimportant 
or where universal wall functions are well established. In 
turbulent boundary layers at low Reynolds numbers, in 
unsteady or in separated flows, or in three-dimensional 
flows, the flow close to the wall must be described. [Ref. 
18] 
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A common algebraic turbulence model divides the flow 
into an inner and outer layer. The inner layer is defined 
by a modified mixing length formula that utilizes some 
damping function. The outer layer includes the wake and 
another damping function. [Ref. 17] 

Other empirical methods currently in use, such as the e n 
method, predict transition. The e n method is a stability 
method, based on linear stability theory. It assumes 
transition begins when a small disturbance is introduced at 
or below a critical Reynolds number. The transition is 
amplified by e 9 . This method allows greater generality of 
the flow, however the formulation still relies on empirical 
terms. [Ref. 19] 

The accuracy of turbulence models are limited by the 
accuracy of the empirical constants. Caution must be taken 
when using a model under different conditions, i.e., a 
different flight regime or a radically different airfoil. 
The turbulence models mentioned above can be fairly simple; 
a more complex model is still not generally usable because 
of the computation costs involved and the uncertainty of the 
constants . 

The influence of turbulence and the transition from 
laminar to turbulent flow on the airfoil need to be 
understood and accurately modeled for a good description of 
the flow over the airfoil to be detailed. As experimental 
methods continue to improve and as computational methods 
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utilize the data, improvement in the detail of the flow 
field and in flow prediction will follow. 



1. Cebeci-Smith Turbulence Models 

The turbulence model used in the Cebeci Interactive 
Boundary Layer (IBL) is a simple algebraic eddy viscosity 
expression. Simple algebraic models seem to adequately 
predict turbulent flow for wall boundary layer flows in 
which the Reynolds shear stress and frequency do not change 
rapidly. However, if the rate of change of shear stress or 
the frequency is large, turbulence models are not currently 
satisfactory. [Ref. 2] 

The Cebeci IBL code utilizes the algebraic 
eddy-viscosity formulation of Cebeci and Smith [Ref. 20]. 
The turbulent eddy viscosity, for wall boundary flows is 
defined by two separate formulas; one for the inner region, 
based on the Van Driest approach, and the other for the 
outer region, based on a velocity defect approach: 




for 0 £ y < y c 




(3.20) 



where : 
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and 
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Y = 



1 + 5.5 (y/S) 2 ' 
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The continuity of the eddy viscosity defines y c ; the 
expression for the inner region is used outward from the 
wall until it agrees with the outer region, which is then 
used. Ytr i s an intermittency factor which allows for a 
transition region when progressing from laminar to turbulent 
flow. It is given by 
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(3.21) 



where the transitional Re number. Re = (u e x/v tr ) . The 

*tr 

empirical constant Gyt r is dependent on the Reynolds number 
of the flow. High Reynolds numbers flows indicate Gyt r = 
1200, lower Reynolds number flows seem to be better modeled 
by lower values of Gy tr . [Ref. 12] 

The parameter a in the outer region is given by 
a = 0.0163/F 2 * 5 where F is the ratio of the normal stress 
turbulent energy to the shear stress turbulent energy 
evaluated at the point of maximum shear stress. This can be 
expressed as 



_ | (u' 2 - v' 2 )3u/3x ) (3 22) 

| - 5^ 3u/3y )<-u'V) max 

The ratio of the time-averaged quantities are 
assumed to be a function of R«p = T v / ( - n ' v ' ) max which can be 
represented by 
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(3.23) 



Thus, the expression for a is 



0.0168 



[Ref. 2] 



a 



[1 - 8 0u/3x)/0u/3y)] 2 * 5 



(3.24) 

Note that the value of Y has the effect of reducing 



the eddy viscosity away from the airfoil surface. This 
turbulence model does not take into account the wake region, 
nor is it validated for separated flow. 

2 . Baldwin-Lomax Turbulence Model 

The Sankar Navier-Stokes (N-S) code also uses an 
algebraic eddy viscosity model, the Balwin-Lomax Turbulence 
Model [Ref. 21]. It is based on the Cebeci-Smith two layer 
model [Ref. 2 2] used in the Cebeci IBL code and may be 
expressed as 



2 




for 0 < y < y 
— — c 




(3.25) 



where 



A = 26. 
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The inner region is the same mixing length formula 
of the Cebeci-Smith model, simplified. No intermittancy 
factor is included (flow is calculated as wholly turbulent) , 
A is a constant rather than being dependent on viscosity and 
velocity gradients, and the velocity profile (du/dy) is 
replaced with the product of vorticity and density. 
[Ref. 13] 

The outer region is based on the wake function, 
F wake and the Klebanoff intermittancy factor, F Kleb(Y) • 
[Ref. 23] 



F «ake = W- 25 W “di/W ' 

W F > “ [1 + 5 ’ 5 1 1 



(3.26) 



(3.27) 



The quantities y max and F max are the maximum values obtained 
from the function 



F(y) = y | a,| {1 - exp (-y + /A + )} (3.28) 

which is a form of the mixing length formula used in the 
inner region. 

u dif the difference between the maximum and 

minimum velocity of the velocity profile. F wa k e is similar 
to the y of the Cebeci-Smith turbulence model and thus also 
reduces eddy viscosity away from the airfoil surface. 
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This model has been used in separated flows and in 
the wake, however its validity is not assured in these regions. 
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IV. THE BOUNDARY LAYER EQUATION 



Prandtl clarified the influence of viscosity in high 
Reynolds flows by simplifying the Navier-Stokes equations to 
yield approximate solutions. He divided the air flow over a 
body into two regions: 

1) The region near the surface where viscous forces 
dominate. 

2) The rest of the flow where inertia forces dominate; 
this region may be considered frictionless and 
potential. 

Consider a 2-D incompressible flow over a body. Most of 
the flow is moving at free stream velocity. However, at the 
surface the velocity is zero, increasing to free stream at 
some distance from the surface as shown in Figure 4.1. In 
this first region, called the boundary layer, the velocity 
gradient normal to the wall, 3u/3y, is very large, as is the 
shearing stress. 



x 




(4.1) 



The two regions are not distinct, but are usually divided at 
the streamline where the velocity reaches 99% of the free 
stream velocity. 

Simplification of the Navier-Stokes equations will be 
accomplished by doing an order of magnitude analysis of each 
term. The following assumptions are made: 
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Figure 4.1 Velocity Profile 
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